#include "Matrix.h"
#include "SLE.cc"


class EigenMatrix 
{
	
public:

	//EigenMatrix():SquareMatrix(){}
	//EigenMatrix(const char* fileName):SquareMatrix(fileName){}

	static double powerMethod(const SquareMatrix& A);
	static double SIPowerMethod(const SquareMatrix& B, double alpha);



};



double EigenMatrix::powerMethod(const SquareMatrix& A)
{
	Matrix x(A.numCols(),1,1);
	Matrix y;
	double c;
	cout.precision(10);
	
	for(long i=0; i< 100; ++i)
	{
		y= A*x;
		c=(y.transepose())[0].largestCoord();
		x=y * (1/c);
		
		
		cout <<"c= " <<c <<"  y = (";
		
		long till=y.numRows();	 
		for(long j=0; j<till;++j)
		{
			cout <<y[j][0] <<"  ";
		}
		
		cout <<"\b\b)\n\n";
	}
		
	return c;

}


double EigenMatrix::SIPowerMethod(const SquareMatrix& B, double alpha)
{
    SquareMatrix A(B);
	
	double c,prevC;
	SquareMatrix I(A.numRows(),0);
	
	long till=A.numRows();
	for(long i=0;i<till;++i)
	{
		I[i][i]=1;
	}
	
	Matrix x(A.numRows(),1,1);
	
	
	A=A-(I*alpha);
	
	SLE sle(A,x);
	
	prevC=0;
	
	while(1)
	{
	    Matrix y=sle.gaussianElimination();
		c= (y.transepose())[0].largestCoord();
		x=y*(1/c);
		sle=SLE(A,x);
		
		cout <<"c= " <<c <<"  y = (";
		
		till=y.numRows();	 
		for(long j=0; j<till;++j)
		{
			cout <<y[j][0] <<"  ";
		}
		
		cout <<"\b\b)\n\n";
		
		if(Vector::nearlyEqual(prevC,c))
		{
		
		    break;
		}
		else
		{
		    prevC=c;
		}
		
	}
		
	return c;

}



/* ##############  End of EigenMatrix methods  ############################# */


int main(int argc, char* argv[])
{
	SquareMatrix A("A6");
	double alpha;
	int garbage;
	
	//EigenMatrix::powerMethod(A);
	for(long i=0;i<A.numRows();++i)
	{
	
	    //find alpha using gerchgorin
		double sum=0;
		for(long j=0;j<A.numCols();++j)
		{
		    if(i==j) continue;
			
			sum=sum + abs(A[i][j]);
		}
		
		cout <<"sum " <<sum <<"   A[i][i]  " <<A[i][i] <<"\n";
		cout <<(abs(A[i][i])-sum) <<"     "  <<(abs(A[i][i])+sum) <<"\n";
		alpha = ((abs(A[i][i])-sum) + (abs(A[i][i])+sum))/2;
		
		double mu =EigenMatrix::SIPowerMethod(A,alpha);
		
		cout <<"Eigenvalue of A for alpha= " <<alpha <<"    is :  "   <<1/mu+alpha <<"\n\n";
		cin >>garbage;
	}

	return 0;

}

